function [Rinc,Rfar,Rlef,FIS,ABS,LK] = resp_kernel_mod(in)
% function [Rinc,Rfar,Rlef,FIS,ABS,LK] = resp_kernel_mod(in)
%
% This returns response functions for a reflector element of one of four
% types: 'biblis', 'zion', 'lra', and 'iaea'.  The two group data from each
% was used to generate responses on 100x100 fine mesh per coarse mesh for
% coarse mesh widths between 3 and 24.
%
% The fits include the size current responses:
%     r11g11, r11g12, r11g22, r12g11, r12g12, r12g22, 
% and the three flux responses
%     phi1g1, phi2g1, phi2g2
%
% All fits are rational with the form
%      response(x) = (p1*x^2 + p2*x + p3)/(x^2 + q1*x + q2)
% and the "fit" tables below have the p and q coefficients for the
% responses in the order listed above.  All fits are very good, with
% R-square = 1.0000 to that precision.
%

Delta = in.xcm(2); 

% biblis
if in.reflector == 0
    
	fit = [ 4.499320530876534e-01     1.076995906790784e+01    -6.208309955216892e+01     2.775482206009754e+01     1.336048170467139e+02   % r11g11
            8.194935099121659e-02     8.142739369160132e-02     1.758770212977565e-01     6.126521878619740e+00     3.490084447346774e+01   % r11g12
            5.629247048266561e-01    -1.339642776540536e+00     1.033093806972309e+00    -7.662326515870207e-01     4.295519412196048e-01   % r11g22
            1.001116555335259e+00    -6.749551039391257e+01     1.279913253088467e+03     7.036545024275460e+02     2.532189531819894e+03   % r12g11
           -3.575493199141279e-03     1.127087447788035e-01     1.252515697531180e-01    -3.381324159563597e+00     4.390947640948107e+01   % r12g12
            2.863137598063410e-03    -1.054759745145268e-01     9.560688947309683e-01    -1.789865615523734e+00     9.005910208432701e+00   % r12g22
            9.525412260068751e-02     7.574064099540371e+00     9.739706865926962e+02     1.991043479450145e+01     9.745357649966929e+02   % p1g1
            2.944488975864534e-02     3.861696917269743e+00    -2.297357776653335e+00     5.236204960578603e-01     1.809773979170768e+02   % p2g1
            2.086462020167209e-02     4.931543308204099e+00     4.289655720521887e+01     5.151708654596831e+00     4.462019903583195e+01]; % p2g2
    R1  = 0.0257622;
    A2  = 0.0715960;
    S12 = 0.0231060;
    
% zion        
elseif in.reflector == 1

	fit = [ 3.557003510951864e-01     2.821156997495499e+01    -1.654847468633183e+02     6.541568444286310e+01     3.401163563425456e+02   % r11g11
            2.140580011822162e-01     2.002482713152605e+00    -2.545541512372851e-01     1.097440798627819e+01     2.571215043553948e+02   % r11g12
            6.058708171560200e-01     3.903138334280130e+02    -3.302672291244962e+02     4.678116341884962e+02     1.076357706465575e+03   % r11g22
            4.526927995518790e-01    -3.047358097377089e+01     5.685632381539262e+02     2.749584686890406e+02     1.147702624273630e+03   % r12g11
           -2.952813616485316e-02     1.066292164689157e+00     8.789911110657859e-01     9.113112411269784e-01     2.300609023108013e+02   % r12g12
            2.256322213925202e+00    -1.504188127317865e+02     2.619258620420118e+03     7.314165705391309e+03     6.533035208163316e+03   % r12g22
            1.288170130489380e-01     4.464256095812315e+00     9.724155686466651e+02     1.883819102225569e+01     9.723998150722979e+02   % p1g1
            3.521729759668158e-01     9.092294143135677e+00    -5.298601848035859e+00    -1.253832296370208e+01     3.794074310685139e+02   % p1g2
            1.604190519322070e-01     3.797836157917877e+00     6.574454471472778e+02     6.920175637866044e+00     6.575969986575503e+02]; % p2g2
    R1  = 0.0295000;
    A2  = 0.0094900;
    S12 = 0.0290300;
% lmw with wrong absorption        
elseif in.reflector == 2
    
	fit = [ 3.255575066765932e-01     2.629576462931440e+01    -1.794825407278231e+02     6.538397823407232e+01     3.661881856865655e+02   % r11g11
            1.991658611892725e-01     1.821152494336344e+00    -1.861984915128856e-01     1.100916830175653e+01     2.509059270361406e+02   % r11g12
            7.706776725209414e-01     1.014851870991747e+02    -7.675873249086234e+01     1.224895110182391e+02     2.623598930868905e+02   % r11g22
            5.826953604840697e-01    -4.007161646281256e+01     7.679015505378665e+02     3.282551534903021e+02     1.547262951817721e+03   % r12g11
           -2.789228966399638e-02     1.012663250460893e+00     7.383782778652370e-01     9.518106724837057e-01     2.211018277678619e+02   % r12g12
            1.259186193895988e+00    -8.210617515211092e+01     1.385371923971967e+03     4.294770538146121e+03     3.454765777954661e+03   % r12g22
            1.285401225246425e-01     4.457740502870401e+00     1.059309808072654e+03     2.051109389815174e+01     1.059293232352996e+03   % p1g1
            3.620199938010021e-01     8.758730547306323e+00    -5.420653951329641e+00    -1.238529895342466e+01     3.686069941840628e+02   % p2g1
            1.555874570926301e-01     3.956303807746599e+00     5.974344585325262e+02     6.749178786390196e+00     5.976562167417494e+02]; % p2g2
    R1  = 0.0302575;
    A2  = 0.0094900;
    S12 = 0.0275969;        
% iaea
else 

	fit = [ 2.537895442920614e-01     7.818738927070152e+00    -8.987239793294324e+01     3.197137853821088e+01     1.815815439088266e+02   % r11g11
            2.652533481194478e-01     2.505377912513875e+00    -1.314012008857986e-01     1.115051834531081e+01     2.494700069052951e+02   % r11g12
            6.994128830885405e-01     1.908155306139513e+02    -1.691881809669867e+02     2.318651698522592e+02     5.387936590495605e+02   % r11g22
            1.920900590349994e-01    -1.305592372358567e+01     2.473663200215944e+02     8.388323258170631e+01     4.995101391078879e+02   % r12g11
           -3.955761197079458e-02     1.426387269713798e+00     9.261528978449024e-01     1.114265209646290e+00     2.042900161490921e+02   % r12g12
            1.225120824138580e+00    -8.330215240022159e+01     1.458593658872336e+03     3.987749695000423e+03     3.431696507623098e+03   % r12g22
            1.186003219252542e-01     4.141001225191800e+00     9.560047922580659e+02     2.323219086613685e+01     9.560669774145848e+02   % p1g1
            4.060619632243836e-01     1.223010008504900e+01    -6.858619875727738e+00    -1.207165285762800e+01     3.763941218718032e+02   % p2g1
            1.583065906995580e-01     3.850811237918328e+00     6.443258287229513e+02     7.071608098840870e+00     6.444814636031144e+02]; % p2g2
    R1  = 0.04;
    A2  = 0.01;
    S12 = 0.04;    
end

        % Responses from fits.
        R11g11 = rational(fit(1,:), Delta);
        R11g12 = rational(fit(2,:), Delta);
        R11g22 = rational(fit(3,:), Delta);
        R12g11 = rational(fit(4,:), Delta);
        R12g12 = rational(fit(5,:), Delta);
        R12g22 = rational(fit(6,:), Delta);    
        phi11  = rational(fit(7,:), Delta);
        phi12  = rational(fit(8,:), Delta);
        phi22  = rational(fit(9,:), Delta); 
        
        % Transverse terms by balance.  
        R13g11  = 0.5*(1 - R11g11 - R12g11 - Delta*phi11*R1                   );
        R13g12  = 0.5*(  - R11g12 - R12g12 + Delta*phi11*S12 - Delta*phi12*A2 );
        R13g22  = 0.5*(1 - R11g22 - R12g22                   - Delta*phi22*A2 );
        
        % Absorption
        ABS(1) = Delta^2*(phi11*(R1-S12) + phi12*A2);
        ABS(2) = Delta^2*(phi22*A2);
        
        % Fission -- none!
        FIS = 0*ABS; 
        
        % Leakage - incident group 1
        LK(1,1) =  Delta*(R11g11 + R11g12 - 1); %incident
        LK(2,1) =  Delta*(R12g11 + R12g12    ); %far
        LK(3,1) =  Delta*(R13g11 + R13g12    ); %to right of incident
        LK(4,1) =  Delta*(R13g11 + R13g12    ); %to left of incident 
        % Leakage - incident group  2
        LK(1,2) =  Delta*(R11g22          - 1); %incident
        LK(2,2) =  Delta*(R12g22 + 0*R12g22    ); %far
        LK(3,2) =  Delta*(R13g22 + 0*R13g22    ); %to right of incident
        LK(4,2) =  Delta*(R13g22 + 0*R13g22    ); %to left of incident 
        
        Rinc(1,1) = R11g11;
        Rinc(2,1) = R11g12;
        Rinc(2,2) = R11g22;
        Rfar(1,1) = R12g11;
        Rfar(2,1) = R12g12;
        Rfar(2,2) = R12g22; 
        Rlef(1,1) = R13g11;
        Rlef(2,1) = R13g12;
        Rlef(1,2) = 0.0;    % no upscatter    
        Rlef(2,2) = R13g22;
        
end

function y = rational(p, x)
y = (p(1)*x^2 + p(2)*x + p(3))/(x^2 + p(4)*x + p(5));
end


